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I. INTRODUCTION 



Green function methods l2| have been very suc- 
cesful in the description of various properties of many- 
electron systems, ranging from atoms and molecules to 
solids [H, [j] ■ Within the Green function approach, these 
properties are completely determined by the self-energy 
operator E, which incorporates all the effects of exchange 
and correlation in a many-particle system [l[ . One of the 
most widely used approximations to the self-energy is the 
GW approximation (GWA) [|. In the GWA, the self- 
energy operator has the simple form £ = —GW, where 
G is the Green function that describes the propagation 
of particles and holes in the system, and W is the dy- 
namically screened interaction. This quantity describes 
how the bare interaction v between electrons is modified 
due to the presence of the other electrons and appears 
as a renormalized interaction in terms of Feynman dia- 
grams. In extended systems the screened interaction is 
much weaker than the bare interaction, and therefore it 
is much more natural to expand the self-energy in terms 
of the screened interaction than in terms of the bare in- 
teraction. The lowest order in this expansion [5[ is the 
GWA. 

Calculations within the GWA are usually done in two 
steps. First, a density functional theory (DFT) [f| calcu- 
lation is performed and the DFT orbitals and eigenvalues 
are used to construct a first guess Go , for the Green func- 
tion and a first guess Wo, for the screened interaction. 
In a second step, the self-energy £ = — G0W0 is con- 
structed and the Dyson equation is solved for the Green 



function. In principle, this new Green function should 
be used to calculate a new self-energy and this process 
should be iterated to self-consistency However, one 
usually stops after the first iteration. The corresponding 
approximation for the Green function is known as the 
Go Wo approximation and has become one of the most 
accurate methods for the calculation of spectral proper- 
ties and band gaps of solids [!, H[. One reason for not 
going beyond the first iteration of the Go Wo method is 
the large computational cost involved. There are fur- 
ther indications that a full self-consistent solution would 
worsen the spectral properties as a consequence of a can- 
cellation between dressing of Green functions and vertex 
corrections Q. This was investigated for the electron 
gas [8j and the Hubbard model [9(. However, this prob- 
lem has not been investigated in detail for real systems 
mainly due to the computational cost involved. 
The Go Wo approximation has, however, two unsatisfac- 
tory aspects. The first aspect is related to the satisfaction 
of conservation laws. Baym [l(| has shown that the self- 
energy expressions that can be obtained as a functional 
derivative of a functional $[G] of the Green function, 
i.e. £ = 5Q/5G, have the important property that they 
lead to conserving many-body approximations. These 
approximations obey basic conservation laws, like the 
ones for particle number, momentum, angular momen- 
tum and energy. The GWA is one of these conserving 
schemes [III [III, [l3j] . However, the ^-derivable approx- 
imations are only conserving when the Dyson equation 
for the Green function is solved fully self-consistently. 
A lack of full self-consistency will generally result in a 
violation of the conservation laws. For this reason the 
use of conserving approximations, such as GW, is crucial 
in obtaining a correct description of transport phenom- 
ena within a nonequilibrium Green function approach 
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EE Oil E3, [3 ■ Since it is one of our research goals 
to study quantum transport, it will be necessary to con- 
sider the fully self-consistent GW (SC-GW) approxima- 
tion @, m m m, m m si. 

A second unsatisfactory aspect of nonsclf-consistent 
schemes, such as GqWo, is that the values of the ob- 
servables depend on the way they are calculated. For 
instance, the total energy can be calculated in different 
ways from the Green function and the self-energy: using 
the Galitskii-Migdal formula [25[ , a coupling constant in- 
tegration [lH , a Luttinger-Ward expression [H, [H, [2?| 
|28[ or various other expressions. For nonsclf-consistent 
calculations all these expressions lead to different results 
and therefore to ambiguity in the value of the energy. 
It was, however, demonstrated in the work of Baym [101 ] 
that self-consistent ^-derivable approximations are not 
only conserving but also have the property that all the 
various ways in which the observables are calculated pro- 
vide the same result. This is another motivation for con- 
sidering fully self-consistent many-body schemes. 
We can therefore conclude that self-consistency is im- 
portant to obtain conserving and unambiguous results. 
However, the large computational cost of self-consistent 
schemes makes them unattractive for the calculation of 
the properties of large and extended systems. In order 
to lower the computational effort it is possible to use 
partial self-consistency which may result in a less severe 
violation of conservation laws. One can, for instance, 
keep the screened interaction fixed during iteration of 
the Dyson equation. This leads to a scheme that can 
be shown to still conserve the particle number and that 
has been tested on the electron gas 0, [2j|. Another 
approach in which the self-consistency is constrained is 
the so-called quasi-particle self-consistent GW (QSGW) 
method [3(| HH, H2, [HJ. In this approach a frequency 
independent self-energy of GW-form is constructed and 
used to solve a quasi-particle equation from which the 
Green function and the screened interaction are con- 
structed iteratively. Due to the Hermitian nature of the 
self-energy the method leads to an orthonormal set of 
quasi-particle states and thereby restricts the form of 
the Green function and the screened interaction. This 
method has been succesful in improving the Go Wo band 
gaps and band widths for a large range of solids [32| . 
One could further consider similar other approximations 
within a quasi-particle framework [34j . Such approxi- 
mations have been shown to improve the band struc- 
ture when local density approximation is a poor start- 
ing point. These methods are, however, not ^-derivable 
and are in general not conserving. Extending methods 
based on quasi-particle equations to the time-dependent 
case is not as straightforward as for the SC-GW, GWo 
and Go Wo methods, which are instead based on an equa- 
tion of motion for the Green function. For the same rea- 
son the computational schemes used in this paper (which 
aims at an extension to the time-dependent case) would 
need to be modified in order to do QSGW calculations. 
We therefore did not consider the QSGW method in 



this work. However, we propose another partially self- 
consistent scheme which is computationally cheaper than 
the GWo method. In this approximation the correlation 
part of the self-energy is fixed during the iteration cycle 
while only the Hartree and exchange parts are updated 
sclf-consistently. In this paper we investigate this ap- 
proximation and other GW schemes at different levels 
of self-consistency and test them on atoms and diatomic 
molecules. We also present in more detail the computa- 
tional method behind the self-consistent GW calculations 
that we described briefly in an earlier Letter [35[ . The pa- 
per is divided as follows: In Sec.|TT]we briefly present the 
general formalism and in Sec. lIIII we describe in detail the 
GW approximation at different levels of self-consistency. 
We then present in Sec. IIVI the details of our compu- 
tational procedure. Finally, in Sec. [Vj we will discuss 
the results obtained with the GWA at different levels of 
self-consistency for atoms and some diatomic molecules. 
These systems are well-suited to test the GW at differ- 
ent levels of self-consistency, but we are ultimately in- 
terested in applications in quantum transport theory for 
molecules attached to macroscopic leads. In such appli- 
cations the long range screening effects, as incorporated 
in the GWA, are important. The investigations in this 
paper are a first step in this direction and aim to get 
further insight into various aspects of the GWA that are 
relevant in quantum transport theory. 



II. GENERAL FORMALISM 

We study finite many-particle systems using the Mat- 
subara formalism [l], l3q ] which can easily be extended to 
a nonequilibrium version of the theory [33, [H, HI] . We 
consider a many-body system in thermal equilibrium at 
a temperature T and chemical potential /z, and with the 
Hamiltonian (in second quantization [l|) 

H = J dxft(x)h{r)$(x) + 

+~ J J dxidx2^ t (xi)^ t (x2)w(ri,r 2 )V'(x2)^(xiX.l) 

Here x = (r, a) denotes the space- and spin coordinates. 
The two-body interaction v is taken to be of Coulombic 
form v(ri,r 2 ) = f/|i"i — r 2 1 - We use atomic units h — 
m = e = 1 throughout this paper. The single particle 
part of the Hamiltonian h(r) has the explicit form 

Mr) = -iv 2 + u,(r)- M , (2) 

where w(r) is the external potential and where we ab- 
sorbed the chemical potential /i into h. The equilibrium 
expectation value of an operator O in the grand canonical 
ensemble is then given by 

(6)=Tv{pd}, (3) 



where p 



- P -0H 



/Tre P H is the statistical operator, 



j3 = 1/ksT the inverse temperature and ks is the Boltz- 
mann constant. The trace is taken over all states in Fock 
space The Green function is then defined as 

G(xr 1 ,x'r 2 ) = - 6( Tl - r 2 )(fe(xr 1 )4(x'r 2 )) 

+ 0(t 2 - ri)<^(x'r 2 )ViH(xT 1 )), (4) 

where we define the Heisenberg form of the operators in 
this equation to be Oh — e TH Oe~ rH . Since the Hamilto- 
nian is time-translation invariant, the equilibrium Green 
function only depends on the difference between the time 
coordinates: G(xti,x't 2 ) = G(x, x';ti — r 2 ). The Green 
function satisfies the equation of motion 

-d T -/i,(r)]G(x,x';T) = 
= <5(t)5(x - x') + 

J dn J dxiE[G(|(x,xi;r-ri)G(xi,x';T 1 ), (5) 

where the self-energy £[G](x, x';t) incorporates the 
many-body interactions of the system. The self-energy 
can be approximated with the usual diagrammatic meth- 
ods Since £[G] is a functional of the Green function 
Eq.© must be solved self-consistently. The self-energy 
is usually split into a Hartree part and an exchange- 
correlation part, according to 

E [G] (xi , x 2 ; r) = S(r)5 (xi -x 2 )v H (ri )+£ xc [G] (x x , x 2 ; r) , 

where the Hartree potential is defined as the potential 
due to the electron charge by 



vh(y) — / dx'n(x')u(r, r'), 



where we introduced the electron density 
n(x) = lim G(x, x; —77). 

77— >Q 



(7) 



(8) 



The main task is now to find an approximation for this 
exchange-correlation part £ xc of the self-energy and to 
solve Eq.©. We convert Eq.© to integral form [|(| 

G(xi,x 2 ;t) = G (xi,x 2 ;r) 

+ J dTidr 2 J dx 3 dx4Go(xi,x 3 ;T - n) 

x (E[G](x 3 ,x 4 ;ri - r 2 ) - S(n - t 2 )£ (x 3 , x 4 )) 

x G(x 4 ,x 2 ;r 2 ). (9) 

Here we introduced a static reference self-energy £0 and 
a reference Green function Go which is defined by the 
equation 



-d T -h(r)\ G (x,x';r) = 

= 5(t)<5(x - x') + J rfx 1 E (x,x 1 )Go(x 1 ,x / ;r)-(10) 



GW 
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FIG. 1: The GW self-energy E is the functional derivative 
of a functional 5>[G]. 



In practice we solve first Eq. (|T0|) for Go and then we 
solve Eq.© for G. It is clear from Eq.([5]) that a fully 
self-consistent solution of Eq. © does not depend on the 
reference Green function Go- In this work we choose 
for Eo a Hartree-Fock (HF) or a density functional self- 
energy. In the first case E = vh [Go] + Ea;[Go], consisting 
of Hartree and exchange parts, whereas in the second case 
Eo = #(x — x')uhxc[Go](x), where uhxc(x) is the sum of 
the Hartree and the exchange-correlation potential [6( . 

From the Green function several observables can be 
calculated. To calculate the total energy E = T + V ne + 
Uq + U xc we use the fact that the exchange-correlation 
part U xc of the interaction energy is given by [l|, [2| 



^ J dr J dyL i y dx 2 E xc (xi,x 2 ; -t)G(x 2 ,xi;t). 

(11) 

The kinetic energy T, the nuclear-electron attrac- 
tion energy V ne , and the Hartree energy Uq = 
1/2 J drdr'n(r)v(r, r')n(r') can all be calculated directly 
from the Green function. To calculate the ionization po- 
tentials from the Green function we used the extended 
Koopmans theorem 0, [H, |H, 13, [H[ , a short deriva- 
tion of which is given in Appendix |E] 



III. THE GW APPROXIMATION AT 
DIFFERENT LEVELS OF SELF-CONSISTENCY 

A. Fully self-consistent GW 

Within the GWA the exchange-correlation part of the 
self-energy has the explicit form 0, [4(| |47j 

£ xc (xi,x 2 ;t) = -G(xi,x 2 ;t)FF(xi,x 2 ;t), (12) 

in which IT is a dynamically screened interaction corre- 
sponding to an infinite summation of bubble diagrams 
(see Fig.[T}. From this figure we see that this self-energy 
is given as a functional derivative of a functional ^[G] 
with respect to G and hence represents a conserving ap- 
proximation From the diagrammatic structure we 
see that the screened potential W satisfies the equation 

W(xi,x 2 ; t) = v(ri, t 2 )S(t) + 

+ J dx 3 dxi J dr'i;(ri,r 3 )P(x 3 ,X4;r-T')lT(x4,x 2 ;(i^) 
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where v is the bare Coulomb interaction and P is the 
irreducible polarization 

P(xi,x 2 ;t) = G(xi,x 2 ;t)G(x 2 ,xi;-t). (14) 



The problem is now completely defined. Equations ijlc 
and Eq. fM]) need to be solved self-consistently together 
with Eqs.JH]), © and 



B. The GoWo and GWo approximations 

The GqWq approximation, as mentioned before, is ob- 
tained from a single iteration of the Dyson equation 
Eq.©, starting from a refence Green function Go- For 
this approximation the self-energy is given as £ xc [Go] = 
— GqWq where Wo is calculated by inserting Go into 
Eq. (fl"4]) and solving Eq. (fl3|) with this irreducible polar- 
ization. The Dyson equation ([9]) is then solved with this 
self-energy to obtain an improved Green function G from 
which spectral properties are calculated. In principle one 
should insert this Green function into the self-energy and 
solve the Dyson equation again for a new Green function. 
This procedure should be continued until self-consistency 
is achieved, but this is rarely done in practice for the rea- 
sons mentioned in the introduction. 
We further consider a partially self-consistent scheme in 
which we write the self-energy as £ XC [G, Go] = —GWo, 
where the Green function G is determined fully self- 
consistently by repeated solution of the Dyson equa- 
tion and where W is calculated from G in the same 
way as for the Go Wo approximation. This reduces the 
computational cost considerably as it avoids the self- 
consistent calculation of the screened interaction W. The 
corresponding approximation is known as the GWo ap- 
proximation [29l . l48| . This approximation was shown 
to be number conserving by Holm and von Barth [491 ] 
for the case of homogeneous systems. More precisely 
they derived that the GWo approximation satisfies the 
Hugenholtz-van Hove theorem [501 ] for the homogeneous 
electron gas. However, one can readily derive the num- 
ber conserving property for the inhomogeneous and time- 
dependent case. This requires nonequilibrium Green 
functions in the proof, but this extension is straightfor- 
ward [51| . If we regard Wo as a given potential (albeit 
nonlocal in space and time), it is clear that £ = <5$/<5G 
for $[G, W ] = -l/2trGGW , where the trace denotes 



integration over space-time variables. Since this $ is in- 
variant under gauge transformations (the phases cancel 
at each vertex of $), we can follow the proof of Baym [l(| 
and derive that GWo is particle conserving. However, for 
time-dependent and inhomogeneous systems Wo is not 
invariant under spatial and time-translations, unlike the 
bare interaction v that usually appears in the functional 
$ [G] . Therefore the GWo approximation will not be mo- 
mentum or energy conserving. 



C. The GWf c approximation 

The most time-consuming part of the GWo calculation 
is the evaluation of the correlation part of the self-energy 
which is nonlocal in time. We therefore propose another 
partial self-consistent scheme in which we only evaluate 
the time-local Hartree and exchange parts of the self- 
energy in a self-consistent manner. We therefore split 
the self-energy as follows 



S[G,Go]-S ffi? [G] + E c [G 



(15) 



The first term in this equation represents the Hartree- 
Fock part of the self-energy 



V Ht [G] = UH [G] + E X [G] 



(16) 



which consists of a Hartree part and an exchange part 
£ X [G] = — Gv. The last term in Eq.([15|) represents the 
correlation part of the self-energy and has the explict 
form 



S C [G ] =-G (Wo-v), 



(17) 



where Wo is calculated from Go in the same way as for 
the GoWo approximation. The approximation for the 
self-energy of Eq. fTB")) will be denoted as the GWf c ap- 
proximation (where fc stands for fixed correlation). This 
approximation is not conserving but, as we will see later, 
nevertheless produces observables in very close agree- 
ment with those obtained from a fully SC-GW calcu- 
lation. 



IV. COMPUTATIONAL METHOD 

A. Numerical solution of the Dyson equation 

In the following, we will describe the computational 
methods that we employed for calculating the Green 
function and the screened interaction W. We consider 
the case of spin-unpolarized systems where the Green 
function has the form 



G(x,x';r) =<WG(r,r';r) 



(18) 



The calculations are carried out using a set of basis func- 
tions such that the spin-independent part of the Green 
function is expressed as 



G(r, 



^G^T^r^V). 



(19) 



The basis functions (f>i are represented as linear combi- 



nations of Slater functions ipi(r) 



which are centered on the different nuclei and are charac- 
terized by quantum numbers (n^, li,mi) and an exponent 
Xi. In these expressions and Y, mi (ft) are the usual spheri- 
cal harmonics. The molecular orbitals 4>i and eigenvalues 
ti are obtained from a Hartree-Fock or DFT Kohn-Sham 
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calculation in this basis. The particle number N is de- 
termined by the chemical potential. Since we consider 
closed shell systems we have N/2 doubly occupied HF 
or Kohn-Sham levels et (some of which may be degener- 
ate) . We therefore choose [i such that e$ = — fi < for 
i < N/2 and a > for i > N/2. In the zero-temperature 
limit (we used j3 — 100) the observables are insensitive 
to the value of fx, provided ejv/2 < M < e N/2+i- The 
reference Green function Go corresponding to the Hamil- 
tonian /i +E (either HF or DFT) is diagonal in the basis 
{4>{\ i.e. in matrix form we have Gijo( r ) — $ijGifi{ T ), 
where 

G i>0 (r) = 6{T)(n(ei) - l) e - £ > T + 0(-r)n( ei )e- e * T , (20) 

and n(ej) — (e^ ej + 1) _1 is the Fermi-Dirac distribution. 
The Dyson equation of Eq. ([9]) in basis representation has 
the form 

G(r) = G (t) + 

dr> f dr"Go(r-r')S c [G,G ](r'-r")G(T"()21) 
'o Jo 

where we denote 



S c [G,G ](r) = S[G](r) - <5(t)£ [G 



(22) 



and where all quantities are matrices. Since in the limit 
t — > — , G yields the density matrix, it is convenient 
to solve the Dyson equation for negative r- values. We 
therefore rewrite Eq. pTj) as 

Gij(r) = SijGifi(r)+ 

[ dn [ d7aG <l0 (r-n)S§ b (r 1 -7a)G? i y(7>0?3) 

with t € [—(3, 0] where we changed variables t\ = r' — /3, 
t-2 = t" — (3, and used Gq(t) — — Go(r + /3) with the same 
relation for G [4fJ. We now discretize Eq. using a 
trapezoidal rule on a time grid (t^ = 0, t^- 1 ' . . . , t^" 1 - 1 = 
—(3). Since the Green functions behave exponentially 
near the endpoints of the imaginary time interval [—0, 0], 
we used a uniform power-mesh [20]. We briefly describe 
this mesh in Appendix [B] The discretized version of 
Eq. ([23]) attains the form 



E 

k.q 



A T (?) 



G w (r<«>), (24) 



where we defined as 



Z ik (r^\r^) = [ dTG lfi (r^-r)j: c lk (T-T^). (25) 

The time steps are positive, where At^ 9 ' = r^ 9-1 )— r^ q+1 ^ 
except at the endpoints where At^ = — and 



Ar (m) = r (m-i)_ r (m) i For a fixed j, Eq. ([24]) represents 
a set of linear equations of the form 



lX3) 



(26) 



where 

A Qi,Qi 



A 



(ip)(kq) — &ik$pq 



A T (<3) 



(j) (i) 

and the vectors Xq 2 , are defined to be 

The self-energy S c of Ea. ([22]) has the form 

S^.(r) = S c ,„[G](r) + <5(r) [£^(0")] - E° ] , (27) 

where Y, HF is the Hartree-Fock part of the self-energy 
defined in Eq. ([T6]) and S C [G] the remaining correlation 
part. The convolution integral (I25p can therefore be sim- 
plified to 

Z a (rW r«) = Gil0 (rW - rW) [S« F [G(0-)] - S° fc ] 
+ / drG il0 (T^ - r)S c . lfc (r - r^)(28) 

When we specify the explicit form of £ c , the solution of 
the Dyson equation is reduced to a calculation of Eq. ([2"5]) 
together with the linear system of equations (|26[) . What 
remains to be discussed is the calculation of the self- 
energy itself. This is discussed in the next section. 



B. Numerical calculation of the screened potential: 
The product basis technique 

To calculate the self-energy we need to solve the equa- 
tion for the screened interaction. The screened interac- 
tion has a singular time-local part representing the bare 
interaction v. It is therefore convenient to subtract v 
from W and to treat its contribution to the self-energy 
explicitly (this is simply the exchange part of the self- 
energy). From the remaining time nonlocal part of W, 

given by W(t x ,t 2 ;t) = W{v 1 ,y 2 ;t) ~ 5(t)v(ti, r 2 ), we 
can calculate the correlation part of the self-energy 



S c (ri,r 2 ;r) = -G(r 1 ,r 2 ;r)W(r 1 ,T 2 ;T). 



(29) 



After this quantity has been calculated it can then simply 
be added to the Hartree-Fock part of the self-energy to 
obtain the full self-energy E[G]. The time- nonlocal part 
W of the screened interaction satisfies the equation 

W(ti,t 2 ]t) = / e&-3«fr 4 «(ri,r3)P(r 3l r4;T)u(r4,r2) + 



dr' / tfr 3 dr 4 t;(r 1 ,rs)P(r3,r4;T - t')W {r^r^) 
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where 



P(r 1 ,r 2 ;r) = 2G(r 1) r 2 ;r)G(r 2! r 1 ;-T). (31) 



The factor of 2 in this expression results from spin- 
integrations in the equation of W using the form of the 
Green function of Eq. (fl8|) . We now insert into Eq. (f3"Tj) 
the basis set expansion for the Green function of Eq. (fT9|) . 
to obtain 

ijkl 

(32) 

where Pijki = 2Gij(r)Gki(—T). By defining the two- 
electron integrals 



W pqrs (T) = J dT 1 dr 2 (f>l(rx)4>*(T2)W(rx,T 2 ;T)(f) r (T2)(f) s (Ti) 

v P qrs = /"dr 1 dr20*(r 1 )0*(r 2 )'u(r 1 ,r2)^ r (r2)^(r 1 ) 
we transform Eq. (|30j) into the equation 

Wp qrs (r) = ^2v p li s P ijk l(T)v jqrk + 
ijkl 

E/ dr' v p i is P ijk i(T ~ T')Wjg r k(r'). (33) 
ijkl Ja 

If we use the multi-indices Qi = (ps), Q 2 = (rq), Q3 = 
(il) and Q4 = (Jk), then we can write this equation in a 
more convenient form as 

Wq 1 q 2 (t)= E v QiQ 3 p Q 3 qA t ) v QaQ2 + 
Q3Q4 

E [ dr'^QsPQsQAr-r'WQ^T'), (34) 
Q3Q4 Jo 

where we defined vu.kj = Vijki and similarly for W and 
Pu,jk = Pijki- We have now obtained an equation which 
we can solve with the same algorithm we used for the 
Dyson equation. 

Note that in this case we effectively use a product ba- 
sis f q (r) — <fri(r)(f>*(r), where q — (ij) is a multi-index. 
This product basis is nonorthogonal and its size is in gen- 
eral much larger than we need in practice due to linear 
dependencies. We thus follow a technique developed by 
Aryasetiawan and Gunnarsson [52j, which allows to re- 
duce significantly the size of the product basis {f q (r)} 
and the computational cost. 

The overlap matrix S for the set of orbitals f g (r) 



On 



(f q \U), 



is diagonalized by a unitary matrix U 

y ] Uq qi (fgi I fq 2 ) Uq 2 qi — UqSqq' , 



(35) 



(36) 



where the eigenvalues a q are positive since S is a positive 
definite matrix. We now define a new set of orthonormal 
orbitals o„ as 



^fE^ww, 



(37) 



with (g q \g q >) = S qq >. Our strategy is use the orbitals g q 
as a new basis and discard the functions that correspond 
to a q < e (we used e = 1CP 6 ). This leads to a much re- 
duced basis as compared to the set of all functions f q . As 
described in Ref. I52L this corresponds to discarding func- 
tions that are nearly linearly dependent and contribute 
little in the expansion. The quantities S, W and P will 
be represented in this new basis using 



(38) 



For the irreducible polarization we then find from 
Eq. {32]) that 



F(r 1) r 2 ;T) = En g '(T)/ g (ri)/ 3 *(r 2 ) 

qq' 



E 



J2 U L P lAr)Uq'q 2 
qq' 



9192 

where 

Pqq,{T)=2G ij {T)G kl {-T). 

With q = (il) and q' — (jk) we have 

P(ri,r 2 ;r) = E ^naflgi ( v i)9*q 2 0*2 ), 



where 



P qiq2 = [y/*U*P(r)Uy/Z\ 



(40) 



(41) 



(42) 



and yfa is the diagonal matrix (\fa) pq — Spq^/Oy. To 
calculate the screened potential we now insert Eq. (j4"T]) 
into Eq. (|30p and readily obtain the matrix product 



W qq ,(T)= VP(T)V 



vP(t -t')W(t') 



qq' 



where we defined the matrices 

W qq ,= J ^T^ng^WiTurr^g^Tn) 

and 



(43) 



(44) 



(45) 



9192 



It is important to note that in Eq. (|4"Tj) and Eq. (|4"3")) the 
summation only runs over the indices q for which a q > e. 
We see from Eq. (l4"2"j) that terms with a q < e contribute 
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little to the total sum. This leads to a considerable reduc- 
tion of the number of matrix elements for v, P and W. 
Finally the correlation part of the self-energy of Eq. (|29|) 
is given by 



i(r) 



kl pq 
•3 



x J ^r 2 ^(r 2 )«^(r 2 )«7*(r 2 ) 
— E Gki(r)Zikji, 



/ ^ra^rOEcfa.rajT^ra) 



(46) 



where 



GW 

gw 

— GW ic 

-- G„w' 



Zik,jl = ^ TjVpUik,pW pq (T)Ul u ^/Oq 



3 4 
R (a.u.) 



FIG. 2: Particle number for H2 at different interatomic dis- 
uy\ tances within the SC-GW, GWq, GWf c and G0W0 approxi- 
mations. 



We can summarize our procedure as follows: in the first 
step the overlap matrix S qq ' of Eq. (f35|) is obtained and di- 
agonalized. Further, using and Eq. (f3"T)) and (|4"S")l the two- 
electron integrals in the new basis Vpq 8X6 constructed for 
p and q such that o~ p , a q > e. Subsequently, for the same 
values of p and q the matrix P pq (r) is constructed from 
Eq.(|l2]) and W pq (r) is solved from Ea.(|33]). In the last 
step, the matrix (|47|) is obtained and the self-energy is 
calculated from Eq. (|4"6")l and further used in the solution 
of the Dyson equation. 



V. RESULTS 

The various GW schemes described in section Mil are 
applied to a set of atoms and diatomic molecules using 
the computational method of section IIVI Details on the 
basis sets are provided in Ref. (53j |. In general we found 
that, in single processor calculations, the computational 
cost of the GWi c method is comparable to that of the 
GqWq method, and roughly twice as fast as the GWq 
method. The the fully self-consistent GW calculations 
were the most time-consuming. 

Particle number conservation. We start by investigat- 
ing the number conservation property of the different 
GW schemes. In Fig. [5] we display the particle num- 
ber obtained from the trace of the Green function for 
the case of the hydrogen molecule H 2 for different sep- 
arations of the nuclei. We display results for the case 
of SC-GW, GW , GW ic and G W , in which the refer- 
ence Green function Go is obtained from a Hartree-Fock 
calculation. We see that the SC-GW and GWq schemes 
yield an integer particle number of N — 2 for all internu- 
clear separations. This is a consequence of the number 
conserving property of both approximations. This can 
be seen as follows. If we would adiabatically switch-on 
the two-particle interactions from zero to full coupling 
strength within a conserving scheme then the particle 



number would be conserved during the switching. This 
is because the conserving property is independent of the 
strenght of the interaction and follows from the structure 
of the ^-functional only. Therefore the particle number 
of the final correlated state will be the same as the parti- 
cle number of the initially noninteracting system. Hence 
conserving schemes always yield integer particle number 
for finite systems at zero temperature. For the case of 
the hydrogen molecule this is N — 2 for all bond dis- 
tances. For the case of GqWq we see that the particle 
number conservation is violated as the particle number 
deviates from N = 2 for all bond distances, the largest 
deviations occuring for the larger bond distances. For 
the larger separations left-right correlation [54[ in the hy- 
drogen molecule, not incorporated in the Hartree-Fock 
part of the self-energy, become increasingly important. 
This puts more demands on the quality of the correla- 
tion part of the self-energy and consequently nonconser- 
vation of the particle number becomes more apparent 
at longer bond distances. Although the violation seems 
small (about 0.01 electron at R — 4.5) it should be em- 
phasized that a change in particle number of 0.05 can 
give large changes in the spectral features and conduc- 
tive properties for molecules attached to leads. A clear 
example of this is presented in the work of Thygesen [l6| . 
For the GW{ C (See sec. IIII CP we also observe a violation 
of the number conservation law with increasing error for 
larger internuclear separations. The error with respect 
to GqWq is however reduced by a factor of 3 at R = 5.5 
as a consequence of a partial inclusion of self-consistency. 

Ground state energies. For the various GW schemes of 
section Hill we calculated the total energies of some atoms 
and diatomic molecules from Eq. (|ll[) . The reference 
Green function Go for the nonsclf-consistent schemes was 
obtained from a Hartree-Fock calculation. In Table U 
we show the results. From comparison with benchmark 
configuration interaction (CI) results we see that the to- 
tal energies of atoms and molecules calculated within all 
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schemes are not very accurate. However, as we will see 
later, energy differences are much better produced. We 
can nevertheless make a number of useful observations 
from the total energies. We first note that all approx- 
imations produce a total energy that is lower than the 
benchmark CI result, with the GoWo generally produc- 
ing the lowest and thereby the worst values. Both the 
GWq and the GWf c methods yield total energies in ex- 
cellent agreement with SC-GW results, where for most 
systems the difference is 10~ 3 Hartree or less. This means 
that both the GW and the GWf c methods can be used 
to make an accurate prediction for the SC-GW energy 
at a much lower computational cost than the fully self- 
consistent calculation. 

Binding curve. The calculation of binding curves is a 
good test for the quality of total energy calculations. In 
Fig. [3] we display the binding curve of the H 2 molecule 
for the various GW schemes together with benchmark 
CI results. The reference Green function Go was taken 
from a Hartree- Fock calculation. We further checked that 
using a Go obtained from an LDA calculation only influ- 
ences the results slightly. For the values of the energies 
around the bond minimum we see the same trend that 
we observed before: all GW schemes lead to a total en- 
ergy that is lower than the benchmark CI results with 
Go Wo being the lowest. The total energies of the par- 
tially self-consistent schemes GW and GWf c are very 
close to the fully self-consistent GW results for all bond 
distances. Although all GW schemes considerably im- 
prove the bonding curve obtained from an uncorrelated 
Hartree- Fock calculation it is clear that all these schemes 
deviate considerably from the CI results in the infinite 
atomic separation limit. To cure this feature one either 
has to do a spin-polarized calculation or go beyond the 
GW approximation and include vertex diagrams in the 
diagrammatic expansion for the self-energy. The shape 
of the binding curve around the bond minimum is well 
reproduced by the SC-GW, GWo and GW/ C schemes, 
implying that these methods may be used to obtain ac- 
curate vibrational frequencies. Since the shape of the 
bonding curve is only determined by total energy differ- 
ences, this already indicates that these approximations 
may perform better in obtaining the energy differences 
than in obtaining total energies. 

Two- electron removal energies. To test the perfor- 
mance of the various GW schemes in obtaining energy 
differences, we investigated the two-electron removal en- 
ergies of the beryllium and magnesium atom. Since these 
atoms and their doubly ionized counterparts are closed 
shell they were suitable test systems. Moreover, the 
beryllium atom is a well-known case for which electron 
correlations play an important role due to strong mix- 
ing of the 2s and 2p states in a configuration expansion. 
In table [H] we display the two-electron removal ener- 
gies for various GW schemes as well as for the Hartree- 
Fock approximation. The reference Green function Go 
is again obtained from a Hartree-Fock calculation. The 
self-consistent and partially self-consistent GW schemes 
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FIG. 3: The total energy of the H2 molecule, as a function 
of the interatomic distance, calculated from the GW approx- 
imation at various levels of self-consistency and CI 

yield results within 0.1 cV from the experimental values 
and considerable improve the HF values that differ with 
more than 1 eV from experiment. The Go Wo approxi- 
mation does not improve at all on the HF approximation 
and gives considerably worse results than the other GW 
schemes. We further see that both the GWo an d the 
GWfc approximations give removal energies that are in 
excellent agreement with the fully self-consistent GW re- 
sults. 

Ionization Potentials. In Table HTT1 we show the ioniza- 
tion potentials obtained with the various GW methods 
for a number of atoms and diatomic molecules. These 
ionization potentials were obtained using the extended 
Koopmans theorem, as explained in Appendix [AJ For 
Go Wo the results shown in the first column were obtained 
by using a reference Green function Go from a local den- 
sity functional (LDA) calculation using the parametriza- 
tion of the exchange-correlation functional due to Vosko 
et al. 59]. In all other cases we used a reference Green 
function from a Hartree-Fock calculation. We see that 
the ionization potentials of fully self-consistent GW agree 
well with the experimental values, the main exceptions 
being the H2 molecule and the Be atom, which show a 
deviation of respectively 0.8 and 0.5 eV. The other par- 
tially self-consistent approaches GWo and GWf c yield re- 
sults that are very close to the fully self-consistent re- 
sults. The GoWo approximation based on the LDA ref- 
erence Green function performs a bit worse than the self- 
consistent GW scheme. For He and LiH there is an error 
of about 1 eV and for Ne and H2 an error of about 0.5 
eV. Performing a GoWo calculation based on a HF refer- 
ence Go instead improves the results for several systems 
but worsens the agreement for H2 which is 1.1 eV in er- 
ror. The dependence on the reference Green function Go 
within the GoWo method is clearly unsatisfactory. The 
partially self-consistent approximations suffer much less 
from this problem. For those schemes we found that 
changing the reference Green function from a HF one to 
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TABLE I: Total energies (in Hartrees) calculated from the GW approximation at various levels of self-consistency compared 
to CI values. 



System 


E " "[GhfJ 


E [GhfJ 


E ic [GhfJ 


^SC 


CI 


He 


-2.9354 


-2.9271 


-2.9277 


-2.9278 


-2.9037 1 


Be 


-14.7405 


-14.6882 


-14.7032 


-14.7024 


-14.6674 1 


Be 2 + 


-13.6929 


-13.6886 


-13.6887 


-13.6885 


-13.6556 1 


Ne 


-129.0885 


-129.0517 


-129.0506 


-129.0499 


-128.9376 1 


Mg 


-200.2924 


-200.1759 


-200.1775 


-200.1762 


-200.053 1 


Mg 2+ 


-199.3785 


-199.3451 


-199.3454 


-199.3457 


-199.2204 1 


H 2 


-1.1985 


-1.1889 


-1.1891 


-1.1887 


-1.133 2 


LiH 


-8.1113 


-8.0999 


-8.0997 


-8.0995 


-8.040 3 



From Ref. [H. 2 From Ref. |Hj]. 3 From Ref. [13 . 



TABLE II: Two-electron removal energies En -2 — Em (in eV) 
calculated from the Hartree-Fock and from the GW approx- 
imation at various levels of self-consistency, compared to the 
experimental values. 



System 


HF 




AE GW ° 


AE GW fc 




Expt. ] 


Mg - Mg^+ 


21.33 


24.86 


22.61 


22.64 


22.59 


22.68 


Be - Be 2+ 


26.17 


28.50 


27.20 


27.61 


27.59 


27.53 



From Ref. y5S|. 



TABLE III: Ionization potentials (eV) calculated from the 
extended Koopmans theorem from various GW approaches. 



Sys. 


Gi LDA) W 


G^ HF) W 


GW 


GW fc 


SC-GW 


Expt. 1 


He 


23.65 


24.75 


24.59 


24.56 


24.56 


24.59 


Be 


8.88 


9.19 


8.82 


8.81 


8.66 


9.32 


Ne 


21.06 


21.91 


21.90 


21.82 


21.77 


21.56 


Mg 


7.52 


7.69 


7.43 


7.38 


7.28 


7.65 


H 2 


15.92 


16.52 


16.31 


16.22 


16.22 


15.43 


LiH 


6.87 


8.19 


7.71 


7.85 


7.85 


7.9 



scheme proposed by us, yield results in close agreement 
with fully self-consistent GW calculations. We further 
checked the number conservation properties of the var- 
ious schemes. The fully self-consistent GW scheme be- 
ing ^-derivable does satisfy all conservation laws, but 
also the partially self-consistent GWq approximation was 
shown to be number conserving. The nonself-consistent 
GqWq and the partially self-consistent GWf c approxima- 
tions both violate the number conservation laws but, due 
to the partial self-consistency in GWf c , the errors are 
much reduced in this scheme. A major advantage of the 
latter scheme is, however, that it produces results that 
are close to the fully self-consistent GW results at a much 
lower computational cost. It will therefore be very valu- 
able to test this method on solid state systems for which 
self-consistent GW calculations are difficult to perform 
due to the large computational effort. In this way it will 
be possible to get further insight into the performance 
of self-consistent GW for a large class of extended sys- 
tems. Work on application of the fully self-consistent 
GW method to transport phenomena is in progress [l8| . 



'From Ref. 



an LDA one, only slightly changes the results. 



VI. SUMMARY AND CONCLUSIONS 

We investigated the performance of the GW at differ- 
ent levels of self-consistency for the case of atoms and 
diatomic molecules. Our main motivation for studying 
fully self-consistent ^-derivable schemes was that they 
provide unambiguous results for different observables and 
the fact that they satisfy important conservation laws 
that are important in future nonequilibrium applications 
of the theory [l8[ . We adressed the question to what ex- 
tent partially self-consistent schemes can reproduce the 
results of a fully self-consistent GW calculation. We 
found that both the GWq method, as well as the GW{ C 



APPENDIX A: IONIZATION POTENTIALS 
FROM THE EXTENDED KOOPMANS 
THEOREM 

Here we give a brief description on the way we extract 
the ionization energies from the Green function using the 
extended Koopmans theorem [U 51, H3, B El- As 
input, this method only needs the Green function and its 
time derivative at t = 0~ on the imaginary time axis. 
We define an N — 1 particle state 



I* 



N-lr 



> = / dxui(x)^(x)K >, (Al) 



where U,(x) is determined by requiring the functional 



E N - 1 \u i } = 



$ Ar - 1 [u,]|i/|$ 



N-lt 



($jV-lM|$iV-lfUi 



(A2) 



which describes the energy of the N — 1 particle system, 
to be stationary with respect to variations in Uj. This 
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amounts to minimizing the energy of the N — 1 system 
by choosing an optimal value for itj. We find 



dx(^^ t (x')[V'(x),/f]|^) W< (x) = 
(E»-E»- 1 ) [ dx«|^t(x')^(x)|^) Ui (x), (A3) 



where the last term contains the density matrix. This 
quantity is easily obtained from the Green function as 

p(x,x') - «|^(x'r)^(xr)K) = lim G(x,x',-77) 

(A4) 

i.e. jo(x, x') = G(x,x';0 ) or = Gy(0 ) in molecu- 
lar orbital basis [40]. Also the expectation value under 
the integral on the righthand side of Eq. (|A3|) , is easily 
obtained from the Green function 

-9 T G(x,x';r)Uo- - «|^V) $(x), A] |<) 

= A(x,x'). (A5) 

In this derivation we used a zero-temperature formula- 
tion but making a connection to the finite temperature 
formalism is straightforward. When we take into account 
that, in the finite temperature formalism, we included the 
chemical potential in the one-body part of the Hamilto- 
nian (see Eq.([2]), then from l|A3[) and (|A5|) we obtain the 
eigenvalue equation 



J dx A(x,x')uj(x) 



N-l 



fi) / dx p(x, x')uj(x), (A6) 



where p and A are calculated according to Eq. (IA4| lA5j) . 
A similar equation for the electron affinities can simi- 
larly be derived starting from an N + 1-state. Since both 
matrices p and A are easily evaluated from the Green 
function, Eq. (|A6j) provides an easy way to extract re- 
moval energies from knowledge of the Green function on 
the imaginary time axis. 



For completeness we mention that the extended Koop- 
mans method also provides a simple way to extract quasi- 
particle or Dyson orbitals [45[ and to construct the Green 
function on the real frequency axis. The Dyson orbitals 
are given by 



Mx) = <$f- 1 |^(x)K) 



dx' <(x')(* Ar |V t (x')V'(x)K) - 
= J dx'p(x,x'K*(x'). (A7) 

In terms of these orbitals and the extended Koopmans 
eigenvalues the hole-part of the Green function is then 
given on the real frequency axis as 



G(x,x';w) = ^2 ■ 



/n(x)/*(x') 



Similar derivations can be carried out for the affinities 
and the corresponding Dyson orbitals from which the 
particle-part of the Green function can be constructed 
on the real axis. 



APPENDIX B: THE UNIFORM POWER MESH 

The uniform power mesh (UPM) [2(| is a one- 
dimensional grid on an interval [0, 0\ which becomes more 
dense at the endpoints. Therefore, it is well-suited to de- 
scribe the Green function on the imaginary time axis, 
since it behaves exponentially around r = and r = ±/3 
[20l |40T ] . The UPM is defined by two integers u and p and 
the length of the interval (3. The procedure to construct 
it is simple: we consider the 2{p — 1) intervals [Q,(3j] and 
\p-Pj,0\ for j = 1, . . . ,p- 1 with fa = I3/2J, anc i divide 
each of these intervals in 2u subintervals of equal lenght. 
The endpoints of all these intervals define our grid which 
has 2pu + 1 grid points. 
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